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ABSTRACT eRe 1 


Manned aircraft that are intended for surveillance or 
to complete a bombing mission will very likely be engaged 
by surface to-air-missiles having guidance systems based on 
infrared (IR) technology. The objective of this study was to 
characterize via Simulation the amount of “cover” that can 
be obtained by dropping from a pre-launched, unmanned 
tactical air launched decoy (TALD) a sequence of pyrophoric 
materials to create an IR cloud, analogous to the 
interference created by microwave chaff, that would protect 
the manned aircraft from the missile. The performance 
analysis is based on a simple reticle based model in which 
the two-dimensional (2D) image is reduced to either a 
composite signal, created by the aircraft, or a composite 
noise, created by the pyrophoric expandable. The analysis 
leads to a computer simulation model producing time and 
space dependent signal-to-noise ratios. It is demonstrated 
that the simulation model can answer questions such as how 
long the materials need to burn, how much intensity 1s 
needed, what wavelength range is most effective, which 
Py cophoric packets should be dropped, and how many. A 
visual model of the time dependent IR pyrophoric cloud has 


also been created. 
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I. INTRODUCTION 


A. BACKGROUND 

The most effective anti-aircraft weapon ever created 
to-date is the infrared (IR) guided missile. Since their 
introduction into operational use in the early 1950s, IR 
missiles have far exceeded the expectations their designers. 
The enabling technology for practical IR missiles came from 
WWII research on detector materials that could be used to 
fashion detectors sensitive in the IR bands. These initial 
materials were sensitive to radiation in the near-IR (1-2 
micron wavelength) band. Over the last 20 years, advanced 
designs have taken advantage of better detector materials, 
moving the bandpass from the near-IR into the mid-IR (3-5 
micron band), where engine plume provides a superior signal 
for all-aspect engagement. The three to five micron band 
also has less atmospheric attenuation and clutter factors. 

The need to protect aircraft from attack by effective, 
easily launched IR homing missiles led to the development of 
pyrotechnic flares. The flare needs to have a jamming-to- 
Signal (J/S) ratio of greater than one to one versus the 
protected aircraft’s signature to lure a missile away from 
the aircraft. Since the observed surface area of the IR 
Signature of the flare is quite small, a pyrotechnic flare 


has to operate at very high temperatures to match or exceed 


the in-band aircraft signature. But high temperatures shift 
the wavelength, and we are forced by the laws of physics to 
increase the size of the flare to provide enough energy at 
the correct wavelength. The flare’s intensity is presumed to 
rapidly rise to several orders of magnitude higher than that 
of the target signature and ignite near the target aircraft, 


inside the field of view (FOV) of missile’s seeker. 


Bs APPROACH 

In this thesis the amount of protection that can be 
obtained by dropping a sequence of pyrophoric materials to 
create an IR chaff cloud will be simulated. A simple 
scenario will be defined in which the packets are dropped 
from a pre-launched unmanned tactical air launched decoy 
(TALD), and the manned aircraft follows, at a later time, 
close to the same path as the decoy, taking advantage of the 
interfering effect created by the pyrophoric clouds. 

The heat-seeking missile 1s assumed to be ain a 
prelaunch status so that missile aero-dynamics are not part 
of this initial study. The imaging capability of the missile 
is assumed to be defined in terms of a definable number of 
pixels within the known field-of-view. Also assumed known is 
the detector sensitivity and detectivity over one or more 
bands of wavelengths. The image of the aircraft is defined 
in terms of a distribution of hot spots. Any interference 


provided by the IR signature of the TALD aircraft are 


ignored. A visual model for expenditure of the IR chaff 
cloud is created. The sequence and timing for the dropping 
of the packets are assumed to be known. Also, both the speed 
of the manned aircraft and the decoy are assumed known. 

The analysis of the problem can answer questions such 
as how long the materials need to burn, how much intensity 
1s needed, what wavelength range and distribution are most 
effective, the rate at which pyrophoric packets should be 
dropped and how many, and lastly, quantitatively what kind 
of signal-to-noise ratio STING, deterioration can be 
extracted. The mathematical software MATLAB will be used to 
Simulate the pyrophoric material, provide image snapshots in 
time of the pyrophoric material, and calculate the resulting 
S/N. 

Some basic IR theory such as Plank’s Law, Wien’s Law, 
and the Stefan-Boltzmann Law, are reviewed in Chapter II. 
Radiometric definitions, atmospheric absorption, and general 
characteristics of the different types of IR detectors are 
also discussed in Chapter II. The design of the pyrophoric 
model is discussed in Chapter III. The design of the 
aircraft plume model is addressed in Chapter IV. The 
integration of the pyrophoric and plume model is discussed 
in Chapter V. The first part of the simulation code and some 
S/N curves are presented in Chapter VI. The second part of 


the simulation code and the pyrophoric and plume images are 


presented iain Chapter VIE; The use of muleaple pyrophonic 
expendables is discussed in Chapter ail. Finaligg 
conclusions from the study of the pyrophoric model and 
future enhancements in modeling are presented in Chapter IX. 
Also, in Appendix A, B and C ais a brief description of the 
symbols and the quantities that are used during the analysis 
or in the simulation code, the MATLAB codes, and the 


mathematical derivations, respectively. 


II. BACKGROUND 


A. INTRODUCTION 

In the year 1800, Sir William Herschel, the royal 
astronomer to the King of England, was conducting an 
experiment with a prism in sunlight. The prism spread the 
sun’s rayS into a spectrum from violet to red. Herschel 
placed a thermometer in the violet color and recorded the 
temperature. He moved the thermometer through the colors 
from blue to red and noticed that the temperature increased 
progressively. He then moved the thermometer beyond the red 
end of the visible region and the temperature continued to 
increase. Thus, he found energy beyond the red; this energy 
has come to be known as infrared (IR). 

A portion of the electromagnetic spectrum is indicated 
in Figure 2.1. All electromagnetic radiation obeys similar 
laws Ot reflection, BeraaGt Lom, diffraction, and 
polarization. The velocity of propagation is the same for 
all. They differ from one another only in wavelength and 
frequency. The portion of the spectrum that includes the 
infrared is depicted in greater detail in the lower part of 
Figure 2.1, where the infrared region is bounded on the 
short-wavelength side by visible light and on the long- 
wavelength side by microwaves. It is convenient to subdivide 


the infrared region into several parts. There are no exact 


designations for the separation of infrared bands. The most 
universally accepted designations today are as follows: 
visible light up to 3 um is designated near infrared (NIR), 
middle infrared (MIR) is from 3 to 6 um, and far infrared 
(FIR) is from 6 to 15 um. Beyond 15 um is the extreme 


Dbmikeared (XGR). 
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Figure 2.1 The Electromagnetic Spectrum. (From Ref [1]) 


B. RADIOMETRIC SYMBOLS 

Ald objects with momzero temperature Willi enuse 
electromagnetic radiation. Some objects prove to be better 
radiators than others. The (hypothetical) best radiator is 
one that obeys a set of classical equations called the 
blackbody radiation equations. Some of terms used to 


describe IR sources and radiation are defined im Table 22s 


Table 2.1 Symbols, Descriptions, and Units 


Radiant Energy transferred 
electromagnetic waves 
Rate of transfer 
































Radiant Radiant flux emitted per 
W-cm™ 
W emittance unit area of a source 
Radiant Radiant flux per 
W-cm 


intensity solid angle 


: 


Spectral 





Radiant flux per Ullal ak Ie, 





solid angle per unit area 






Radiant flux incident per 








unit area 





Radiant emittance per unit | 








radiant wavelength interval at a 






emittance | particular wavelength 








Ratio of radiant emittance. 
of a source to that of a 
E | Emissivity blackbody at the same| Dimensionless 


temperature | 
i 


(Cs BLACKBODY RADIATORS 
The spectral density distribution of the power radiated 
into a hemisphere from a blackbody radiator per unit area of 


the source iS given by Planck’s Law [2]: 


w, (a Hoo ] cy ] 
soy 2 explic/,-)-1 cae x “Yar )- (11-1) 


where 


-W, 1s the spectral radiant emittance of the source with 
units of W-cm~*-um”™ 


-h is Planck’s constant (h=6.625x10~° joules:s) 
-c is the vacuum velocity of light (c=3x10° ms”) 
-i is the wavelength of the radiation (um or nm) 
-k is Boltzman’s constant (k=1.38x10” joules-K") 
-T ais the temperature of the source in K 

-¢,= 2mthe’=3.74x10° Wem *-um' 


we.— Mey k-1 439x100 salen-K 


Plots of the spectral radiant emittance for a typical range 
of temperatures encountered in our environment are shown in 
Figure 2.2. The spectral radiant emittance is the spectral 
power distribution normalized by the area of the source. 
Note that most of the energy is in the infrared portion of 


the spectrum. As the temperature of the source increases, we 


note that the area under the curve increases, and the 


location of the peak shifts toward the shorter wavelengths. 
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Figure 2.2 Spectral Radiant Emittance vs. Wavelength on 
Blackbodies of Various Temperatures. (From Ref [2]) 

The location of the peak of the spectral radiant emittance 
is found by taking the derivative of Planck’s Law with 
respect to wavelength. It can be shown that the result gives 
Wien’s Law, the wavelength of the peak spectral radiant 


emittance 


_ 2.89610" 
ae iT 7 [uum] (II - 2) 


Knowing the temperature of the source, we can calculate the 


location of the peak emittance and, as shown in Figure 2.2, 
higher temperatures imply lower j.. 

The value of the peak spectral radiant emittance is found by 
Suisstttuting, (il=2) "ante: 4011-1 )\eto cbhotaam 


4 ie S 
WAC et ae (II -3) 


CA 


where 

b =1.286x10°° W-cm™- um" -K>. 
The total emittance W is the total power emitted into a 
hemisphere at all wavelengths per unit area of the source. 


Ie is found by integrating Planck’s Law over all 


wavelengths, 
w=|W,da (11-4) 
0 


The result of this integration is 


IIe 


= = IT-5 
15c*h? ( 


where 
Oo =5.67x10°°W -cm? °K 


which represents the Stefan-Boltzmann Law and shows that the 
total power per unit area of the source increases as the 


fourth power of the temperature. 
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We frequently pass electromagnetic radiation through a 


spectral filter (like the atmosphere) with a spectral 
response of the form S(A). The emittance that is passed by 


sucn a £fiiter is 
Ris 
W(A,,As)= fs(aW,aa (1-6) 
a 


Figure 2.3 is an illustration of one such graph where the 


vertical axis is used to compute the fraction of the total 
emittance that lies between a wavelength of O mm and some 


upper wavelength divided by the total emittance over all 


wavelengths. 
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Figure 2.3 W(0,A) Relative to Total Emittance. (Ref [2]) 
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D. RADIATION FROM OBJECTS 
In order to calculate the radiation emitted from a 
source, we model the source as either a uniformly radiating 
spherical source(e.g, the sun or a light bulb) or as a 
planar source(e.g, the exhaust port of a jet engine). 
ae Uniform Spherical Source 
For a spherical source with a known temperature(in K), 
an emitting area A, and an emissivity(usually assumed to be 
less watlan ) One; 1-627 aS Gravybody sob ject [2], we can 
calculate the radiant emittance of the source as 
W = eoT* (II -7) 
Then we can find the total power from the source as 
P=W-A (1 -8) 
and the radiant intensity by realizing that it 1s emitting 


uniformly into a sphere of 4” steradians: 


ye 


oS IIT-9 
re (11-9) 
2 Planar Source 
For a planar source, we cannot assume that the 


radiation 1s symmetrical SO the power 1s radiated 
nonuniformly from such a source. A standard procedure is to 
model the source as a Lambertian emitter, where both power 
and radiant intensity will vary with the observation angle. 


For such a source it can be shown that the radiance is [2] 


y -W _ sor" (11-10) 
7 Tt 7 HG 
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E. TRANSMISSION OF IR THROUGH THE EARTH’S ATMOSPHERE 

Most infrared systems muSt view their targets through 
the earth’s atmosphere. Before it reaches the infrared 
sensor, the radiant flux from the target is selectively 
absorbed by several of the atmospheric gases, and scattered 
away from the line of sight by small particles suspended in 
the atmosphere. 

When the particles are small compared with the 


wavelength of the radiation, the process is known as 
Rayleigh scattering and exhibits a A“ dependence. For larger 


particles, the Scattering is independent of wavelength. 
Scattering by gas molecules in the atmosphere is, therefore, 
negligibly small for wavelengths longer than 2 um. Smoke and 
light mist particles are also usually small with respect to 
infrared wavelengths, and infrared radiation can, therefore, 
penetrate further through smoke and mists than visible 
radiation. However, rain, fog particles, and aerosols are 
larger and, consequently, scatter infrared and visible 
radiation to a similar degree. 

In the infrared portion of the spectrum, the absorption 
process poses a far more serious problem than does the 
scattering process. The spectral transmittance measured over 
a 6000 ft horizontal path at sea level is shown in Figure 
2.4. The molecule responsible for each absorption band, 


either water vapor, carbon dioxide, or ozone, is shown in 
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the lower pare [OL ~ehneu: Our Ce IMSPeCrICieOrmeamynnC 224 
reveals the presence of atmospheric windows, i.e, regions of 
reduced atmospheric attenuation. 

IR detection systems are designed to operate in these 
windows. Combinations of detectors and spectral bandpass 
filters are selected to define the operating region to 
conform a window to maximize performance and minimize 


background contributions. 
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Figure 2.4 Transmittance of Atmosphere over one Nautical 
Mile (NM) Sea Level Path (Infrared Region). (From Ref [3]) 
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F. IR DETECTORS 

The detector is the heart of every IR system because it 
converts scene radiation into some other measurable form; 
this can be an electrical current, or a change in some 
physical property of the detector. 

In general, an infrared detector in a particular set of 
operating conditions is characterized by two performance 
measures: the responsivity R and the specific detectivity 
D*. The responsivity is the gain of the detector expressed 
in volts of output signal per watt of input signal. The 
specific detectivity is the detector output signal-to-noise 
ratio for one watt of input signal, normalized to a unit 
sensitive detector area and a unit electrical bandwidth. The 
different types of IR detectors may be divided into two 
broad classes, namely thermal detectors and photon, or 
quantum, detectors. 


5 Thermal Detectors 


Because thermal detectors have been used since 
Herschel’s discovery of the infrared portion of the 
spectrum, it 1S appropriate to consider them first. They are 
distinguished as a class by the observation that the heating 
effect of the incident radiation causes a change in some 
physical property of the detector. 

Since most thermal detectors do not require cooling, 


they have found almost universal acceptance in certain field 


TRS 


applications in which Jt as <impwacticalmiteowpeovicde such 
cooling. Because (theoretically) they respond equally to all 
wavelengths, thermal detectors are often used 2 
radiometers. The time constant of a thermal detector is 
usually a few milliseconds or longer, so that they are 


rarely used in search systems or in any other application in 
which high data rates are required. 


2 Photon Detectors 


Most photon detectors have a detectivity that is one or 
two orders of magnitude greater than that of thermal 
detectors. This higher detectivity does not come for free, 
however, Since many photon detectors will not function 
unless they are cooled to cryogenic temperatures. Because of 
the direct interaction between the incident photons and the 
electrons of the detector material, the response time of 
photon detectors is very short; most have time constants of 
a few microseconds rather than the few milliseconds typical 
of thermal detectors. 

Finally, the spectral response of photon detectors, 
unlike that of thermal detectors, varies with wavelength. 
Curves of D* versus wavelength are shown in Figure 2.5 for 


different types of photon detectors. 
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Figure 2.5 D* vs. xX for Representative Detectors (From Ref 


[2]) 
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III. PYROPHORIC MODEL 


A. GENERAL DESIGN REQUIREMENTS 

Infrared decoys are used in modern times to protect 
combat systems from IR tracking threats. The most common 
examples of IR decoys are IR flares. IR decoys are used to 
protect aircraft from heat-seeking missiles. Some key 
requirements for the design specifications of an expendable 
decoy are discussed in the following. 

= i Peak Intensity 

Peak intensity is normally the most TMmoOr taint 
requirement. IR decoys must radiate with sufficient 
intensity and at least exceed the intended target’s radiant 
intensity in the band of interest. This can be accomplished 
by controlling the decoy temperature or by the use of 
selective emitters that have a higher emissivity in the band 
of interest. A common design problem is illustrated in the 
spectral intensity versus wavelength plot of Figure 3.1. The 
relative spectral distribution between a decoy (a small hot 
source) and an aircraft target (a large, relatively cool 
source) is shown. From the viewpoint of the decoy designer, 
it is advantageous if a significantly higher percentage of 
the optical signal at the detector, within the ‘ decoy band’ 
is due to the decoy rather than the target. Finally, the 


peak intensity is the primary driver for the decoy weight, 


19 


volume, and cost that must be smaller, lighter, and 


cheaper,respectively, than the protected target. 
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Figure 3.1 Typical Decoy and Target Spectra (From Ref [4]) 


2% Flare Intensity Rise Time 


A decoy must persist long enough to ensure that no 
possibility of target reacquisition remains. Consequently, 
it must maintain a credible signature until the original 
target is no longer in the threat field-of-view. If this is 
not the case, it is necessary to deploy a second flare. 
Additionally, the decoy must achieve an effective intensity 
quickly enough to capture the seeker before leaving the 


threat field-of-view. The diameter of the threat field-of- 
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view at the time of decoy deployment is usually less than 
200 m. This means that effective operational intensity 
levels must be achieved in a fraction of a second [Ref 4]. 
The exact value of the time to rise to peak radiant 
intensity depends of the chemical composition and packing. 


Exact values are difficult to find in the open literature. 


B. DEVELOPMENT OF THE PYROPHORIC MODEL 

The pyrophoric model must simulate changes in both the 
time and space domains. The modeling procedure for 
calculating the in-band radiant intensity is summarized in 


Beegure 3.2. 





Specify 
Pyrophoric 
Inputs 






Calculate emittance Calculate pyrophoric 
over all wavelengths total radiant intensity 


Calculate emittance 
between the two Calculate in-band radiant 


wavelengths intensity of pyrophoric 
material 


Calculate pyrophoric 
area at max radiant 
intensity 





Figure 3.2 Calculation of In-Band Radiant Intensity. 
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dig Time Domain Variations 
The emittance over all wavelengths for any hot object 
follows Planks radiation equation is given from the Stefan- 


Boltzmann Law 


4 
W(A,T ) op = EOT py (1 -1) 
where the Stefan Boltzmann radiation constant is [Ref 2] 
o =5.67X10°°W -cm? °K™ 
€ is the emissivity of the pyrophoric and T,, a specified 


temperature. The approximation of a graybody model is used 
to approximate the in-band radiant emittance for these 
expendables [Ref 4]. 

The M-file plankpf.m recorded for reference in Appendix 
B is used to calculate for a particular temperature the 
integrated emittance over a band of wavelengths. The 3-to 5- 
micron band is very common for a detector because this band 


has less atmospheric attenuation. The emittance W(A,,A,),, 


that lies inside the detector’s band is obtained using the 


M-file program cited above. The fraction of the emitted 


power that lies between two wavelengths A, and A, is 


w(A,,4), : 
~— (A,,,) "Ay an, (11 - 2) 


The pyrophoric source is modeled as a uniformly 


radiating spherical source with radius r,,, at the instant 
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LOGe vt iehwcmmaxi MUM OGeULrsS Inene sadiant jantensaty t.. The 
associated area 1s 


(III - 3) 


gh 
Mn c= AE Tj ae 


Pf 


The peak radiant intensity J,, of an isotropic radiator is 


(Ref 2] 


J sux = a) ae (11 - 4) 


The total radiant intensity that lies inside the detector’s 


band is 


NCEE) = Jax “Mpg (I ~5) 


ISi2Onit Ne Ge “ab je ie jel selehe tela jepineicjelolonenel™ leybe 9,7 ielev= alias 
radiant intensity distribution of the pyrophoric burn can be 
predicted. Two operational properties are desirable. 
Ideally, the pyrophoric time function will achieve peak 
intensity quick enough to capture missile seeker attention 
before leaving the field-of-view. Furthermore, it should 
persist long enough to ensure that no possibility of target 
reacquisition remains. The probability density function 
(pdf) of the standard Gamma distribution [Ref 5] has 
sufficient flexibility to model these two basic properties. 


The flow chart that summarizes the development of the 
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temporal variation in radiant intensity distribution of the 


pyrophoric burn is shown in Figure 3.3. 








Set 
Pyrophoric 
Inputs 


Create pyrophoric time function 
array 


Calculate parameter a of 
Gamma distnbution for 
pyrophoric time dependence of 
radiant intensity 


Calculate Gamma function for 
pyrophoric time dependence of 
radiant intensity 








Create the time radiant intensity 
distribution of pyrophoric burn 


Figure SIGS: The Time Radiant Intensity Distribution 
Structure Of fy eoOpne re se lisne 
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OEher inputs Vsuch asythat corresponding to initial time t, 


and final time t, for the simulation to be performed will 
lead to the development of a time-function-array t. 


The Gamma distribution, which is a continuous random 


variable (RV), is defined as 





7 ee ec, 7 
fa. B)= Fay e/® (II - 6) 
t20 
a>0 .99) 0 


Because of interest in modeling the temporal variation in 


radiant intensity, the t parameter is used to define time 
dependence. The standard Gamma distribution [Ref 5] has B=1 


so the pdf of a RV modeled as a Standard Gamma RV is 


ee (111 —7) 


t=0 
a>0OQ 


In Appendix C, a short proof is presented that demonstrates 


that the peak in the distribution satisfies the condition 


t =a-l (111 - 8) 


p 


From this relation the values for the parameter a for 
different values of the pyrophoric peak times are 


predicted. Using the M-file gammaz.m, recorded for reference 


in Appendix B, we calculate the Gamma function [(a) and 


2 


employ the main program to predict the temporal variation in 


the radiant intensity as 


(III -9) 


J{t)= J ax  f(t,a) 


11sec, 


D 


Conditions €t 


the 


equation (III-9) for 


Ot 


AS plot 


1.0m for a pyrophoric time array lasting 


ONC. 


=2000K, 


pt 


at 


8s are shown in Figure 3.4. 
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Ze Space Domain 


The flow chart that is followed to develop the 
mic atmelTi tens) tywecdustributiconmeot GsynOphen:c bunnieane the 


space domain is shown in Figure 3.5. 


Represent in space each 
instant of the pyrophoric 
function using a 
Gaussian distribution 


Calculate the percentage 
of each Gaussian 
distribution which is not 
on the negative axis side 


Normalize each curve and 
set the values of the mean 
and standard deviation 


Create the weighted values 
for each Gaussian 
normalized space 

distribution 


Integrated view model 





Figure 3.5 Normalized Space Distribution Structure 
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Tt is assumed that the spatial distribution for each 
instant of the pyrophoric burn is a Gaussian distribution. 
Generally, a continuous random variable X is said to have a 


Gaussian distribution if the pdf is 


f(siMo)= poe — ee ~coc ces (111-10) 
7 -O 


where W is the mean and o is the standard deviation. One 
Gaussian distribution is shown in Figure 3.6. The Gaussian 
pdf is symmetric about the mean and the standard deviation 


is the distance from the mean to the point at which the 


slope of the curve changes. 


f(x) 





Figure 3.6 Graph of a Gaussian Distribution 
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In Appendix C a proof is presented that demonstrates that 
for a Gaussian random variable Z, the area under the curve 
from zero to infinity is given by equation (C-12), which is 


reproduced here for convenience: 


ate=-4)= Yen“. 75) (a1) 


Powe ew eles fF, Fi ., and Ene jpycphoriteGmeime | lumection 


array t, the following parameters can be obtained: 


UL, = igh (constant pyrophoric velocity) (111 -12) 
P 

Tyax, =U,*t (temporal dependent pyrophoric radius) (11 -13) 

Tran = -MAX (normalized pyrophoric radius) (II -14) 


MAX 


At each instant of the pyrophoric burn, it is assumed that 
the spatial radiant intensity function iS a Gaussian 


distribution with a mean value 


= itp OF (111 -15) 
Tax aes l, 
From (III-11), the area under each curve from zero to 


PHirinicty is 
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ant | 


~ ft 
Or) = V, -erfe( V5 -a[D (III -16) 


which iS positive and less than one. 


For each element of the array oer the value of the 


temporal distribution is given by 





- eet ei 
a) = toe III -17 
fle.a)= 7 (1-17) 
Now in equation (III-10) we use x= fa” and Ve aes ES ObEaum 
the normalized Space distribution in the pyrophoric burn 


as 


_ oF {rm A)/. II -18 
1 pene t ; C)= e 20 ( 


p 21 -O 


The next step is to use 1/Q(t) to normalize the corresponding 


equation (II1I-18). The result, which is a pdf defined over 


positive x has an area under the curve equal to one. 


NORM ? 


This result is multiplied by equation (III-9) to obtain 
Fon, '/ |= fC LY aI Ox Yorg _ 
{ NORM x) F norm a )xJ(t)x( olf) (II 19) 


where the space integration of the space-time distribution 


satisfies the relation 
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This forces the resulting space-time distribution j to have 
a space integrated area equal to J(t). One example is shown 
in Figure 3.7 where successive curves represent time steps 


oer 0.28. 
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Figure 3.7 Pyrophoric Space Distribution 
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The last major step, see the flowchart in Figure 3.5, is to 
create the space domain pyrophoric integrated view-model. 
The general idea is to pick a curve from Figure 3.7 which 
corresponds to one instant of time t, as illustrated in 


Figure 3.8. 
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"NORM 
Figure 3.8 Pyrophoric Space distribution for one instant of 
ete 
As represented on Figure 3.9, the pyrophoric source is 
modeled as spherical. It is assumed that an external remote 
observer will ‘view’ the source and perceive radiation 
integrated along a thin line x. The total radiant intensity 


distribution at a specific instant of time will consist of 


the integration of radiant intensity over each distance Vee, 


where Ir is the distance between the center of the 


+ 


a2 


distribution and some point along the line x. The line x is 


shifted by je] from the distribution center. 


view 





Figure 3.9 The integrated pyrophoric model projection 


Analytically, for each distance from the center of the 
cliches Ve thwalteral | xl, a fixed number of uniformly spaced points 
are selected along the line x from x=0 to x=x,,. The 


distance | that connects each of these points with the 


center of the distribution is given by 


al lz, +|z)° (III — 20) 








l, 
where | x,] the distance along the line x. In the computer 


code, an index k is generated as 
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= ne by somu ) x size(F omy) (111 - 21) 


so the integrated-view radiant intensity for the specific 


lx, 1-6, ONG pecinet is Given by 4G ik) he summation oLecu 


a 


the radiant intensity values for each | r,| give the 


integrated-view radiant intensity 


= x ‘ ° 
Jel »y Se) k index (III — 22) 
which is easily computed. Alternatively, an equivalent but 


more analytical representation is 


: t 
Jy (F,.1)= ee il Fo 
oe a ale (ieee) 


Typical distributions of integrated-view radiant intensities 
for pyrophoric burns lasting eight seconds are shown in 


Figure 3.10. The horizontal axis of Figure 3.10 is specified 


in terms of the index variable of the array r.. It 1s 


observed that the curve with the maximum area corresponds to 


the instant ES. . 
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3.10 Integrated-view Radiant Intensity Distribution 


Curves for the Pyrophoric Model. 


Figure 


33) 


7 





IV. PLUME MODEL 


Ae GENERAL DESIGN REQUIREMENTS 

The plume of the aircraft can dominate the IR 
Signature. This is no doubt due to the considerable radiant 
energy created from the combustion process. Measurements of 
the radiation from military a aircraft’s plume are highly 
elascsificdwmmortumately, it is not toomtdititiveult Yeo gmake 
order-of-magnitude calculations of the radiation by using 
only temperature and dimensional information that can be 
found in the open literature [Ref 4]. 

During afterburning, which results in an increase in 
thrust, the rate at which fuel is consumed increases 
drastically and the plume becomes the dominant source. Also, 
if the aircraft is viewed from the forward hemisphere or 
from an aspect angle at which the tailpipes are not visible, 
the plume is the only source of radiation available. The 
calculation of the exact radiant intensity of the plume is 
extremely difficult since both temperature and emissivity 
vary in a complex manner through out its volume. 

The exhaust temperature contours for a typical turbojet 


engine with and without afterburner are shown in Figure 4.1. 


a7 


260°C ° WITHOUT AFTERBURNER 
he Thrust 15,800 Ib 







— = 
ro] © Sa) 
PEUETTITIty rity | 
ee) 
~~ 
— 
6 
@) 
Ce) 
iy 
2° 
‘@) 
en) 
re) 
© 
© 


Oo oo 


50 100 150 200 250 300 350 400 







Radial distance from engine centerline (ft) 


Distance from face of 
tailpipe (ft) 
5 93°C 66°C 
538°C 
10 427°C 
WITH AFTERBURNER 
Thrust 23,500 |b 

15 

20 


Figure 4.1 Exhaust Temperature Contours for a JT4A Turbojet 
Engine with and without Afterburner (From Ref [1]) 
We can see that when afterburner is turned on, the 
temperature and size of the plume increase appreciably. If 
we integrate over the entire plume in order to estimate its 
radiant intensity, the radiant intensity will be several 
times that of the hot tailpipe. For in-band engineering 
calculations, a plume can be approximated as a graybody with 
an emissivity of 0.9 [Ref 1]. 

The steps that are followed to calculate the radiant 


intensity of the plume model is shown in Figure 4.2. 
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Specify 
Plume Inputs 





Calculate 
emittance over 
all wavelengths 


Calculate 
emittance 
between the two 
wavelengths 


Calculate plume 
area 


Calculate plume 


Calculate plume 
radiant intensity 


total radiant 
intensity 


inside the 
detector’s band 





Figure 4.2 Calculation of In-Band Radiant Intensity. 
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B. DEVELOPMENT OF THE PLUME MODEL 

The assumed geometric model for the plume, shown in 
Figure 4.3, is an ellipse. This is used to model a three- 
dimensional radiating ellipsoidal source. The required 
Paw~ameters eres the included) majormeaxlcwa., Jeane miner vac 


b 


i) and the temperature %) “Cheats 1S Ss essumed ) constam- 


throughout the plume. The area of the ellipse is 


A, =%-a,,-0,, iV =") 





Figure 4.3 Plume’s Radiating Source Area. 


The Stefan-Boltzmann Law iS again used to predict the 
radiant emittance. Here it 1s applied to the case of the 


plume: 
W(A,T) = €0T »' (7V -2) 


where € here is the effective emissivity of the plume. 


Using the M-file plankpl.m i@ag the Darti cu las 


temperature with the specific wavelengths of the detector’s 


band, we can compute the total emittance W(A,,A,) that lies 
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inside the detector’s band. Next, we find the fraction of 


the emittance that lies between wavelengths A, and A,as 


WA, Ay) 
Mas 7 eal 7), (Iv -3) 


The physical model for the plume source is assumed to behave 


as a Lambertian emitter [Ref 2]. Following the standard rule 


for calculating radiance from total emittance [Ref 2], we 
get 
Eo -T;, 
NA ae 
I (IV —4) 


The total radiant intensity is then given by 


Jp, = Np, * Ap, (Iv -5) 


and the plume radiant intensity that lies inside the 


detector’s band is 


Tey cul | ' 
IV -6 


For later reference, the composite image signal associated 


with the plume target is defined 


S=J(A,Ay)p UV -7) 
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V. INTEGRATION OF PYROPHORIC/PLUME MODELS 


The main program, recorded in Appendix B, consists of 
two parts. The general structure of the main program is 


represented in Figure 5.1. 


Specify Inputs 










Generate S/N Ratio with 
variation in parameters 
over range of time 


Generate images at 


one instant 





Figure 5.1 General structure of the main program. 


Once the inputs have been specified, a selection in the 
main program permits execution of either the first part, 
which is the generation of the signal-to-noise ratio (S/N), 
or the second part, which refers to the generation of the 


pyrophoric flare and plume images. 
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The main steps of the procedure which are followed to 


create S/N are shown in Figure 5.2. 


Change 
parameter 





Produced time 
radiant intensity 
GLSErI DUEL On er 
PY Tr ORnOri Gc umm 


Calculate in-band 
radiant intensity 
Of PpVrOpnos1 Cc 


material 


Calculate in-band 
radiant intensity 
of plume 


Generate S/N 


Figure 5.2 S/N Structure. 
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A. CALCULATION OF S/N 

hae rirsee step "Of the flowchart wim Figure S22, the 
computation of the time radiant intensity distribution of 
the pyrophoric burn, was covered in Chapter III. One typical 
distribution is represented in Figure 3.4. 

In the second step of the flowchart, the in-band 
radiant intensity calculation of the pyrophoric material for 
each instant of time t is performed. This is found by 
performing a multiplication of the temporal variation in the 
radiant intensity J(t), given by the equation (III-9), with 


the fraction of the emitted power y,, that lies inside the 


detector’s band (A,4,), given by equation (III-2). 
Next, the in-band radiant intensity of the pyrophoric 


material for each instant of time t is obtained from 


NAb eae = J(t) “Tp (V ~1) 


and 


N=J(A,A,)py re) 


where N is interpreted as a composite noise interference for 


the missile tracking system. 
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The in-band radiant intensity calculation of the plume, 
as given by equation (IV-7) in Chapter IV, is the last 
component needed in the S/N prediction. 

For the calculation of S/N, a situation similar to that 
shown in Figure 5.3 is assumed, where the plume represents 


the signal and the pyrophoric the noise. 


(SIGNAL) 





(NOISE) 


Figure 5.3 S/N Representation. 


Based on the preceding discussion, by taking the ratio of 
the result of equation (IV-7), which represents the 
composite plume image Signal, with the result of equation 
(V-2), which represents the composite pyrophoric image 
Signal, we obtain an ssc iede O/T: 

(A,,A 


J\A,; 2 ye (V 7 3) 
HES Ns 


» 
N 
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After the estimation of S/N, the main program permits a 


change in parameters such as pyrophoric temperature T,,, 
plume temperature T,,, detector’s band wavelengths (i,,A,), 


and a recalculation 


EyEopnomre peak sere Cl, Or (raciucmer 


with the new inputs produces a new value for S/N. 

The second part in the main program refers to the 
generation of the pyrophoric and plume images. The main 
steps of the procedure which are required to create the 


images are shown in Figure 5.4. 
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Calculate in-band 
radiant intensity 
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material 







Create 
normalized space 
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Create integrated 
pyrophoric model 














Create basic 
pyrophoric image 










Calculate in-band 
radiant intensity of 
plume material 







Create basic 
plume 
image 







Normalizations to 
produce radiant 
intensity for the 
pyrophoric and plume 


Generate images 






Figure 5.4 Generate Images Structure. 
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B. GENERATION OF IMAGES 

Ls Pyrophoric Image Construction 

The first four steps of the flowchart in Figure 5.4 
have been already created. Specifically, this includes the 
in-band radiant intensity calculation for each instant of 
time t. This calculation is done in the first part of the 
main program which also deals with the S/N generation. Also 
covered in these four steps is the conversion from the 
normalized space distribution to the integrated pyrophoric 
model discussed in Chapter III. 

To construct the pyrophoric image, it is assumed that a 
Situation similar to Figure 5.5 is applicable, where FOV 
represents the field-of-view of the incoming missile, R is 
the distance at which the missile will detect the pyrophoric 


the number of pixels that have to be used in 


x 


source, and N,, 
order to visualize the pyrophoric image. These quantities 
are some of the inputs which are specified at the beginning 
of the main program. 

From this assumption, the missile’s window dimension 


Foes weaueculated as 


Reg = FOV (rad )x R(m) [m] (V -4) 
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pix 








vy 





2 


Figure 5.5 Pyrophoric Image Construction. 


To specify the distance of the pixels from the center of the 


pyrophoric image, the pythagorean relation is used 


Ne iy (V —5) 


where the i, j indices take values from -(N,,,/2) to (N,,,/2). 
Once N. has been calculated for all pixels, ehie 
corresponding radial distances are predicted from 


pix 
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Next, in the algorithm an index is generated in order to 


assign a value from the integrated-view radiant intensity, 


equation (III-23), to each pixel. This index is determined 
by 
Roy ; 
f= ——— x size(ryon, ) (V -7) 
(r MAX ° max(r NORM ) 


For each value of this index that corresponds to the inside 
of the integrated-view radiant intensity array index, the 


corresponding radiant intensity value is assigned 


C; = J, (é) (V -8) 


otherwise the zero value iS given to the pixel. Finally, a 
normalization is performed. The value of each pixel is 
divided by the summation of all the pixels and multiplied 
with the in-band pyrophoric radiant intensity, equation(V- 


1), (A,,A,)..- The normalization rule is therefore equivalent 


|e) 


ce 
Cy =e x dI(A Ay)», (Vv =8) 


© DC, 





where the time dependence is implicit. 
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Be Plume image construction 

From the flowchart in Figure 5.4, the only necessary 
information to create the basic plume image is’ the 
calculation of in-band radiant intensity of plume material, 
which has been already estimated. Specifically, the in-band 
radiant intensity calculation for each instant of time t is 
given from equation (IV-7). 

To construct the plume image, it is assumed that the 
geometric conditions represented by Figure She: are 


applicable. 
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Figure 5.6 Plume Image Construction. 


Sl 


Note the main difference between Figure 5.5 and Figure 5.6 
is that the plume is represented here as an ellipse instead 


of as a circle. The FOV, R, and N.,. are the same quantities 


that already have been specified for the pyrophoric image 
CONSELrUCETON:. 


To specify the location of the pixels in the plume 


image, one vector r 1s created 


(V-10) 


peat ye. 
where 
- i,j denote indices that take values from -(N,,,/2) to 


(N.,/2) with a specific step 


“ce scenete “Unit vectors along “Ehe ymagor axis ay pane 
minor axis 5, 


- Aa—Red (V -11) 


In order to create an ellipsoidal shape, a value has to be 
assigned to the pixels that belong inside the ellipse’s 


boundary 


(V -12) 








(V -13) 


A 


Taking the preceding into account and making use of 


equations (V-11), (V-12), (V-13), we get 


b a 
Sot | | Pet (Vv -14) 
Riera R rieta 
which provides a convenient rule for defining the interior 


part of the plume in the computer model. 


The numeric values which are assigned to the pixels 
that define the in-band plume radiant intensity J(A,,A,),, are 


taken from plume model and equation (IV-6). After the 
generation of the ellipsoidal shape, a normalization takes 
place. The value of each pixel is divided with the summation 
of all the pixels and multiplied with the in-band plume 


radiant intensity J(A,,A,),,.. This normalization rule is 


equivalent to 





SIERO Wag6), 


FE 


The similarity of equation (V-15) to equation (V-9) for the 


pyrophoric image 1s apparent. 
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VI. PRESENTATION OF S/N CURVES 


As mentioned in Chapter V, the first part of the main 
program, which is recorded in the Appendix B, is primarily 
concerned with the generation of S/N. Once the inputs have 
been specified, the value of S/N can be extracted for each 
instant of time. A temporal variation in S/N during the 
pyrophoric burn time can then be obtained. Next, by changing 


parameters such as pyrophoric temperature T,,, detector band 


Jee Ce BP reg am 


(1,,4,), pyrophoric peak time t,, or radius r 
computes S/N curves for the new inputs. 

A selection of these curves which have been generated 
through the procedure described are presented in this 
chapter. A graph that presents an overall variation of S/N 
Guring the pyrophoric burn time and a chart that records the 
minimum S/N for different parameters are included in each 
figure. 

Specifically, Figure 6.1 is a plot of S/N as a function 
of the pyrophoric function time for different values of 
ie Oriente Mrenperatuee we] sa Resultismawere Generated) Lon 
i OO Ciewmer-l sec, ancdir..—24m As expected, an increaseman 
ao clOniemmuenperatume  T resuilts) -in absdecrease ins /N 


during the burn period of the pyrophoric material. This is 


an expected result because it is known from the Stefan- 


DD 


Boltzmann Law that the emittance increases as the fourth 
power of the temperature. 

Figure 6.2 is a plot of S/N as a function of the 
pyrophoric function time for different values of radius r,, 
corresponding to the instant t, for which a maximum occurs 
in the radiant intensity. Results were generated for 
T,.=2000K, T,,=1000K, and t,=lsec. We observe that increasing 
the radius r,, Significantly decreases S/N between the 
interval and final ~trme Of tHe pyrophoric fEunction. “Pinus 
leads to the conclusion that the radiant intensity output of 
Ene pyrophoric flare increases with r,,.. 

Pagure 6.3 4S a pret oF S/N as "am function of irhe 
pyrophoric burn time for different detector bands. Results 
Were Genereecmuom l=) 00K ame OOK ten = oSECy and aaa ie 
The curve with the lowest values corresponds to the 1-2 
micron band. This band exhibits a large amount of 
atmospheric attenuation and is not a practical window of 
detection for a missile. The next band is the 3-5 micron 
band, which has less atmospheric attenuation. This band, 
along with the 8-12 micron band, exhibits the lowest amount 
of atmospheric attenuation and are practical windows of 
detection for a missile. Making a comparison between these 
two bands, we see from Figure 6.3 that the curve with the 
lowest values corresponds to the 3-5 micron band. The next 


lowest band is the 7-9 micron band. Its curve has lower 
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values of S/N compared with the 8-12 micron band, but again 
this band exhibits a large amount of atmospheric attenuation 
and 1s not a practical window of detection for a missile. 
Figure 6.4 is a plot of S/N as a function of the 
Pipe riortce sence tomucimne Tor Cdirrereneayobues OL tf - Results 
mene Generareq) Lor 1 -2000K, TF j=L000K, and 2 -2m. Form this 
case, the minimum S/N corresponding to edenme were sete ssame- 
The focus here is the length of time for which the 
pyrophoric flare will achieve intensity that can capture the 
missile’s seeker. From this plot, we conclude that the flare 
achieving a peak radiant intensity at a later point in time 
will have an extended windows of effectiveness. This point 
1s demonstrated in Figure 6.4 by comparing the length of 
time S/N is below 0.5. The pyrophoric flare reaching its 
peak at 1.4 seconds has a window of effectiveness of 3.7 
eeeonds, while. themelorer wre = Oo so ponly ettective stor 


2.9 seconds. The former is “better” by approximately 25%. 
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Figure 6.2 S/N for various r,,, 
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Figure 6.3 S/N for various Detector Bandpass 
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VII. PRESENTATION OF IMAGES 


The second part of the main program, which is recorded 
in Appendix B, consists of the generation of the pyrophoric 
and plume images. The procedure which is followed is 
described in Chapter V. 

The image visualization is created by using a definable 
number of pixels. The command “image” in MATLAB creates an 
image graphics object by interpreting each element of the 
pyrophoric or the plume matrix which has been created as an 
index into the figure's colormap. Each element of the above 
matrix specifies the color of a rectilinear patch in the 
image. Along with the command “image”, the command 
“CdataMapping/scaled” is used to scale the values. 

The variation in radiant intensity for the pyrophoric 
image can be described better with two colormaps. The “jet” 
colormap, which ranges from blue to red and passes through 
the colors cyan, yellow, and orange. The “gray” colormap 
returns a linear grayscale colormap. In Figures 7.2 to 7.6, 
a colorbar is used to show the current color scale. 

The second part of the M-file pyrof.m, which is 
recorded in Appendix B, is used to provide image snapshots 
for different instances of time for the pyrophoric model and 


for the plume model which is constant with time. For example 
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see Figure 7.1 which was generated with T,,=1000K. A sequence 
of pyrophoric images on a “jet” colormap for a pyrophoric 
with a burn function of seven seconds, burn temperature of 
i=c000K, detector band 3-5 micron) band, pea tadtaun 
intemesty at t -1.5 seconds, and) radius@ae t equal tomec 
meters are shown in Figures 7.2 through 7.6 for times t=1, 
t=1.5, t=3, t=5, and t=7 seconds, respectively. Looking at 
these plots, we observe that the ring with the maximum 
radiant intensity will be small and close to the center for 
small times. At the end of the burn, for the pyrophoric 
flare the circular periphery of peak radiant intensity has a 
larger size but a lower radiant intensity value. 

One byproduct of the time-space modeling of the 
pyrophoric burn is that the radiant energy profile does not 
remain peaked at the center. This point was previously made 
in Chapter III and qualitatively represented with the time 
dependent graph of Figure 3.8. The images of the pyrophoric 
flare seen in Figures 7.2 to 7.6 demonstrate this property. 
This is a byproduct of the assumption that the peak in the 
ersisbuabele al ghercsgloukicn, aL Ii eleleibuc lee te el ac@ie el Sjolveeuie sl Vela IL 
at u=n, 70. “his 1s “Conststene "wren Ga physica) mode mOrn arm 
expanding material for which the center burns first and then 


an outward radial wave of combustion follows. 
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Figure 7.1 The Plume Image. 
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Figure 7.2 The pyrophoric image at t=1 sec. 
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Figure 7.3 The pyrophoric image at t=t.=1.5 sec. 
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Figure 7.4 The pyrophoric image at t=3 sec. 
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Figure 7.5 The pyrophoric image at t=5 sec. 
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Figure 7.6 The pyrophoric image at t=7 sec. 
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Representation of the images may give the impression that 
setting the plume in the center of the pyrophoric image will 
allow the missile to detect the plume, which is uncovered as 
the radial wave of combustion moves outward. This is not 
true since the reticle based detection system, represented 
in Figure 7.7, converts all the pixels inside the FOV to one 


value. 





Figure 7.7 Reticle based Detection System. 


It should be noted that the S/N predictions in Chapter 
VI are based on this model. In particular, using a reticle, 
the missile electro-optic (EO) process converts  two- 
dimensional (2D) optical image to a one-dimensional (1D) 


electronic signal [Ref 1]. 
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VIII.MULTIPLE PYROPHORIC FLARES 


The main computer program code also simulates the 
dropping of multiple pyrophoric flares, each at different 
instances of time and with generally different parameters. 
Pile Sep areset Choice luce Oo VGophor Cc. temperature le oreadius 
of pyrophoric at the peak radiant intensity, r and the 


MAX, 
ime OF (peal agian Intenienty eat. 

Using this option, a scenario can be generated where a 
pre-launched unmanned tactical air launched decoy (TALD) 
drops a sequence of flares in advance of a manned aircraft. 
The missile will then be distracted by a flare generated IR 
noise curtain. The intent is to provide ‘cover’ to the 
manned aircraft which flies into a safe-zone corridor a few 
seconds later. This is represented in Figure 8.1. 

The settings on the different parameters of the ejected 
pyrophoric can be chosen according to the conclusions 
reached from studying Figures 6.1, 6.2, and 6.4 in Chapter 


Vil. im, goes 1culan, an increase in Chemo GOD nO Lac 


temperature te One an increase vig the radius e 


pf MAX 


Significantly decreases S/N between the initial and final 
time of the pyrophoric function. Also, a flare which 
achieves peak radiant intensity at a later point in time 


will have an extended window of effectiveness. 
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The radiant intensity for four pyrophoric flares, which have 
a function burn time of 5s and are launched at times t=0s, 
t=0.5s, t=ls, and t=1.5s are shown in Figure 8.2. The first 
flare has the shortest time to peak radiant intensity, 
@-0-8s. This value is selected™in order to capture™ the 
missile seeker quickly. Other specifications for the first 
Flare include r,.-2m™and 7 =2200K, whieh were chosen to keep 
S/N low. The second and third flares, which are launched at 
t=0.5s and t=l1s, respectively, have a longer time to reach 
Deak Padtane intenstey, (ures tn sounder temecreare fan 
extended window of effectiveness. For these two flares, 
Get atic 2 OOK ane used. Meserdave been "selected | one. 
Gotoh on ee shee aG, fom Nee race klane Sree WeHere mrss 
still significant radiant energy from the first flare. The 
last tvawe 1s y launched “atist=l.o5 wien larger values Or i. 
and T eseleaciverrereric premlous  eworm i lares, ir =2m, Gauc 
(ace OC hIGhemavyaniecare Useduii -erderm) FO Mmalmeaim a low 
value of S/N by compensating for the reduction in 
effectiveness of the first flare. In addition, t,=1.2s is 
used in order to have an extended window of effectiveness. 
The composite radiant intensity of the four flares 
Versus ene cOtal function our Eime is chown in. Fiocure o.oo. 
The changes in slope are correlated with the activation of 


new flares. 
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The composite S/N versus the total function burn time 
is represented in Figure 8.4. Also, in this plot the changes 
in slope are correlated with the activation of new flares. 

The S/N curves that correspond to each flare are shown 
in Figure 8.5. By comparing Figure 8.4 and Figure 8.5, it is 
obvious there iS an extension of the length of time S/N 
remains below some threshold (e.g., 0.5) with the use of 


multiple flares. This point is quantified in Table 8.1. 





Figure 8.1 Pre-launched Unmanned Tactical Air Launched Decoy 
(TALD) drops a Sequence of Flares in advance of a Manned 
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Figure 8.5 S/N ratio Curves for Each Flare. 


Table 8.1 Time below the 0.5 threshold. 





Time below the 0.5 threshold (sec) 


Flare number 


un) 
7 cA] LO 
d dr] 2 


Composite flare 


Flare 1 
Flare 2 
Flare 3 
Flare 4 





ips 


In order to demonstrate how these results apply to 
tactical decisions that need to be made in the field, two 
examples are presented. The conditions that the manned 
aircraft be given safe cover while being viewed by a 


“stationary missile” can be stated as 


Ate, 2 tage + “4 (VIII -1) 
N ac 


where 


-At..~, 1S the time S/N is depressed below a specified 
Ehnsesnold 


-t,,. 1S the approximate time after the dropping of the first 
flare that the aircraft enters the safe corridor 


-V airerart velocity 


ac 


Equation (VIII-1) can be rearranged as follows: 


= 


l 
|| [ae VIII -2 





For convenience of interpretation, minimum values’ for 
(Wok) are plotted in Figure 8,6 versus At. fOr various 


values in the time delay, t 


aac ° 
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t ] 

step | | 

| step delay [ 0.5 1.0 18 2.0 2.5 3.0] ! 
bain : 


Vac/Rf 





elapsed time S/N depressed 


Figure 8.6 Minimum Velocity Curves for various Values of 
Time Delay Parameter. 


Example 1 


Given V,.=186 m/sec(=0.56 MACH1), R,=62 m, t,,.=1.0s, 
what is the minimum allowed time to keep S/N depressed in 
order to provide “safe” cover? 
Step 1 Calculate V,./R,= 3s” 


Step 2 For delay t,.=1.0s, read from Figure 8.6, 1.5s. 


dac 


From the above graph and reference to Table 8.1, it is 


obvious that any one of the flares studied would be 


U4 


sufficient. However, if the aircraft delay is changed to 3.0 


secomds end V»-62 m/sec, Ehenerme condition oOne=minimum Ata. 


changes to about 4.1 seconds and more than one flare packet 


will need to be dropped. 


Example 2 


Given Ate =2S, R62) m, t..=12 5s, 


(S/N) 
what is the minimum velocity of the aircraft in order to be 


provided “safe” cover? 


Step 1 Ror At 
(V,./R,) Pe 0K 


Bae co wene: je —"obs;, ~cad feom. Higume: 8 jcr 


—~— 
5 — 
min 


SEeCpie2, USeMthemecc TEs mErom stCDelewltl. 2 —o2 om EO Obtaum 
(V =124 m/sec (=37% MACH1) 


ac) MIN 


The point of this example is to demonstrate that a time 
limit in the required reduction of S/N forces a predictable 
lower limit on the velocity of the manned aircraft. 

From Figures 8.4 and 8.6 several significant and 
general observations can be deduced. First, because the 
period of time in which S/N stays below a given threshold is 
finite, as can been seen from Figure 8.4, there 1S an upper 
limit on the delay that can take place before the manned 
aircraft enters the safe corridor. The longest this delay 
can be is the total time interval that S/N is below the 
threshold. Note from Figure 8.6 that each curve predicts the 
minimum required velocity of aircraft which approaches 


infinity as the elapsed time for S/N to be depressed 
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asymptotically approaches the delay. For example, the curve 
with a delay of three seconds approaches infinity from the 
right as the elapsed time S/N depressed approaches three 
seconds. 

These curves predict in general that a greater delay in 
the aircraft entering the corridor requires a greater 
velocity of the manned aircraft in order to guarantee safe 
passage. Conversely, the curves when applied in reverse 
suggest a minimum required time for the depression in S/N. 
In particular, from Figure 8.6 it follows that the greater 
the delay and/or the slower the velocity of the manned 
aircraft, the longer S/N needs to stay depressed below a 


specified threshold. 


qe, 





IX. CONCLUSIONS AND FUTURE ENHANCEMENTS 


The objective of this research was to model the impact 
of dropping pyrophoric type flares which are in the field- 
of-view of a prelaunch missile. It is assumed that these 
flares are launched from an unmanned tactical air launched 
decoy (TALD). The manned aircraft, which needs protective 
cover and will presumably benefit from the flares, will 
follow into the corridor a few seconds after the start of 
the flare sequence. 

In order to reach this objective, a simulation model of 
the time-space radiant intensity distribution for a 
pyrophoric flare expendable has been developed. The 
performance of the pyrophoric flares to create a distraction 
for the missile was characterized in term of S/N where § is 
the radiant power of the plume and N is the radiant power og 
one or more flares. Several conclusions from this study are 
described in the next two paragraphs. 

The radiant power for the pyrophoric flare burn in the 
[3-5] micron band was found to be significantly higher than 
what occurs in the [8-12] micron band. It was also found 
that increasing the pyrophoric temperature Mo he Sihevera Yelejalig\s) 
the radius yr, Significantly decreases S/N between the 
initial and final time of the pyrophoric function. Finally, 


a flare which achieves a peak radiant intensity later in 


Sal 


time has an extended window of effectiveness; i.e., a longer 
window in time for which S/N remains below a specified 
threshold. 

In this thesis the concept of dropping of multiple 
pyrophoric flares to create a “safe” corridor for manned 
aircraft missions was evaluated. In particular, multiple 
flare “packets” lengthen the window over which S/N remains 
below an acceptable threshold. 

Lastly, it was noted that if the manned aircraft is not 
moving with high enough velocity there will be a problem in 
generating “safe” cover. From this concept, a rule relating 
minimum aircraft velocity to both the delay in entering the 
safe corridor and the time S/N is depressed was presented 
and illustrated with several examples. 

A number of simplifying assumptions were applied in the 
study. Future work on the simulation model could focus on 
one or more of these assumptions. For example, to mention a 
few points, atmospheric absorption was neglected in this 
study. Also, the plume was assumed to be a constant 
temperature graybody. A better model for the atmospheric 
absorption would be to consider the cold CO, absorption 
spectrum. For the plume, consideration of the hot CO, 
emission spectrum would be also beneficial (the red spike- 
blue spike effect). Also, an experimental data based model 


for the pyrophoric burn radiation spectrum could be used to 
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refine the simulation model. Finally, the calculations of 
S/N could be based on more complex image processing 
algorithms that are indicative of modern reticle based 


missiles. 
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APPENDIX A. COMPUTER AND ANALYSIS VARIABLES 


Computer Analysis Brief Description 
Variables Variables 


OST TE Se TOS eT ee OE 8 Oe 8 Oe bee ORES GOSS ES BO SE SESS DOMME SESE twas See Sees a NF RS el Et SD SSS ESE Sats DOES TBA eet ee DES SOS TSS HER EST Fe Oe TOSSES 





Missile field-of-view 





Range between missile and target 


PORE be os © ee ee SEE Ee ewe Be SS Oe ee ee eee BES SS ee STE EE See ee SS 


Range inside missile’s field-of- 
view on target’s plane 


Pyrophoric temperature 


ee ese eee See eee eng es peewee 2 oe = |e ewes ee ee 5 5 ee Fe eS em eo ee eee ee ees 





Plume temperature 


Ce ee eee 


Initial wavelength 





Final wavelength 


PP Sees SSeS eS SS TES ee ee ee ee ee ee ee eo owe © ew 2 eee ee Ss te ee ee eB EES ES SSE ESOS SS 


[nm2tral time Of pyLOopnoric 
Ffuneéi on 





Final time of pyrophoric function 


SRO i es ie eee es eee eee ewe es eee ee eee ee ees ce eee ee oe SS ee Se eee ee Se SSeS ee LOSS SSR ERE y DOSS SSS SSS SS Se See Ses Bee ee ee a ee BO ee ee ee ES 


Pyrophoric time at max radiant 
intensity 


Ore. ft eee ee eee ee PS ee ee ee 


Radius of pyrophoric 
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Afl 


Apl A 


pl 


we eee Se Se ee eee eee ewe ees esse ese ss wena s eee es se tester ee esos et ee eee ee ees 


a ae 


eee ee eee errors 


Wftlare 


Sore mes eee ee eee ceo 


ratwtl 





Radius of pyrophoric at max 
radiant intensity 


8 SSS OS 8 OS Se 28S SOS OS = 6 ee ee eS 8 OE Be Oe eee met TOS SO ESS SS 8S ee wee wee eS ee ee eee wee Oe Pee eet eee 


Normalization radius of 
DY rOpaomi 


wee ee ee ee ees ee we ewe ewes Se ee eee tee ee 


Area of pyrophoric at max radiant 
intensity 


OS SSS 28 84 88 ES Fee oe 2 EOE Eee eee ee es ees Cees eee oe oe 222 O89 2S 2S ee Oe 2 8S 2 eee eee eee 


Area of plume 


SO 8 SSF Ce OSS Oe a Se i ESE PS SSeS 


Major axis of plume’s ellipse 
area 


Minor axis of plume’s ellipse 
area 


s eeaieetieniaatiasmatnetinettaehamethednetnatoteteuete titan tietaatietenketaateteatetett bette eet eee He 6 tee ewer ts 


Velocity OL. DyYrOpioric 


Average emissivity 


2 Se eee ose ee we mew Es ee Roe ee ees ee ee Se eee eee 


Stefan-Boltzmann constant 
Total emittance of pyrophoric 
Total emittance of plume 


Sem 88 8 ee 8 ee te et Oe eo i 2 SS OO 8 OE SES See SETS CESS Se Se eee TES 


EM@ecance Of PyropnhoOric in ene 
band of interest 


oe eee we we 8 ee 2 2 ee er a a Oe 8 re ee 8 8 8 8 28 © ee eee eee 


Emittance of plume in the band of 
interest 


SD OS SS EE OS SE 8 8 OS 8S 8 ES SS 8 EE OS 8 OS 8 OE ee ee ee ee eee 


Rractilon sor py roOpnorlG emittance 
in the band of interest 
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Fraction of plume emittance in 
the band of interest 


— ee 


Peak radiant intensity of 
bynegicrts ~~ 
Radiant intensity of pyrophoric 
rune elem 


Oe es we Oe ee © Be 8 oe Pe we 8220S £8 FOOSE S FOOTE Se OT FT © SS eS SS DST SST ET SH SS SESH HS 





Radiant intensity of pyrophoric 

in the band of interest 
Radiant intensity of plume in the 
band of interest 


8 See we wee Bw Oe ee BOS See Bee Se Fe eee ose 


Gamma distribution 


ee re et ee et es ee ESSE wee ee 8 OS Eo es eS we eS SS EE SESS SS 


Parameter (a) of Gamma 
Pi Sea ole oma. 


e882 Boe ee ere we oe ee we emer ee 





Gaussian space distribution 





oon SS SSS £88 OS Se Se ye Se EE SE OT ES OE PRES wee eT SSE we ee ee et ee eee 


Standard deviation of GausSSian 
distribution 


ee a a OSCR eS Se BW ewe ewe Se ee BO Se Ee Pe ees BSS ee eee Fe eS eS SES 


Mean of Gaussian distribution 


Amount under a Gaussian 
Gistribution curve 
Weighted value of a Gaussian 
distribution in normalized space 

| domain 





Integrated view radiant intensity 





SH oe eer eas ewe eee 2 ee he eee meee Sf Oe eee t= 





Number of pixels 


22 2 eee ee Se SS ee = 8S 8 Se ew Sees e868 6 8 8 See Pe See Ee eS Sm a eee ea eee eee 


Pyrophoric image matrix 
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Normalized pyrophoric image 
matrix 


Normalized plume image matrix 


S/N ratio 
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APPENDIX B. MATLAB CODES 

This appendix contains listing of all MATLAB input 
files that were used to get the results posted in Chapters 
VI,VII and VIII. These M-files are: 

1. Pyrof.m, 

2. Plankpf.m 

3. Plankpl.m 

4. Dquadi.m 

5. Dquadl.m, and 
6. Gammaz.m 

Pyrof.m is the main computer code which has two parts. 
The first part generates the S/N curves for one or multiple 
pyrophoric with different parameters and the second part 
generates the pyrophoric and plume images. 

Plankpf.m is a function file which calculates the 
Spectral Radiant Emittance (SRE) for the pyrophoric at 
Pempemcdeane i (Kelvin): 

Plankpl.m is a function file which calculates the 
Spectral Radiant Emittance (SRE) for the plume at 
Eempemaeumet 5 (Kelvin). 

Dquadi.m is a function file which calculates the 
Spectral Radiant Emittance (SRE), for the plume at 


Pom eae Vi Melvin); Simone Ene wccrecrens band: 
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Dquad2.m is a function file which calculates the 
Spectral Radiant Emittance (SRE), for the pyrophoric at 
Gempesatume (Kelvin), (7 incide twerdeteerencs wand. 

Gammaz.m is a function file which calculates the GAMMA 


FUNC LON - 
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Program 1 Pyrof.m 


SESESSEEESE MAIN COMPUTER CODE S8%ESEEEE% 
%$$EtRef. files: 
S plankpft.m, plankpl.m, 
% dquadl.m, dquad2.m, gammaz.m 


format long 


global Tpl 
global lambdal 
global lambda2 
Girooa LOE 
Global sel 
global t2 


$EEE% Inputs B&BS 


program=input(’For the S/N part enter 1, else for the images part enter 
ae iG 

SESESESESESESEESEESESSESE First Part 

ier Ocglham—— | 

Tpl=input(’Set the Temperature of plume (Kelvin): ’); 

Tp=input(’Set the Temperature of Pyrophoric (Kelvin): ’); 
lambdal=input(’Set the value of the lambdal (um): ’); 
lambda2=input(’Set the value of the lambda2 (um): '); 

tl=input(’Set the initial value of the time for the Pyrophoric (sec): 
is 

t2=input(’Set the final value of the time for the Pyrophoric(sec): '’); 
t_peakl=input(’Set the value of the t_peak (sec): ’); 
rmaxp=input(’Set the value of the rmax (m): '’); 

t3=input(’Set the increment step between tl and t2: ’); 

tfl=input(’Set the dropped times of pyroforic: '); 


femtede: tS: EZ: 
pence (1) 
fete 2—1- length (tll); 
ei lame= mila \n 
if length(t_peakl)==1; 
t_peak=t_peakl; 
else 
t_peak=t_peakl(z); 
end 
ar length (hey—— > 
apr=Tp, 
else 
aoe — Loz 
end 
if length (rmaxp) ==1; 
rmax=rmaxp; 
else 
rmax=rmaxp (Z) ; 
end 


; emissivity of Pyrophoric 
7*10%(-12);% Stefan-Boltzmann constant (cm) 


oa 


Wpfb=dquad2 (’plankp£f’,lambdal, lambda2);%Emittance between the two lambda 
%( W cm*(-2) ) 

Wflare=s*Tpf*4; %Stefan-Boltzmann(total emittance over all wavelengths) 
Afl=4*pi* (rmax)*2; % Area Pyrophoric at t_peak 

s1=5.67*10%(-8);% Stefan-Boltzmann constant (m) 

J_max=((e*sl*Tpf£*4)/ (4*pi))*A£f1l; % Total Radiant Intensity(W/sr) 
ratwtl=(Wpfb/Wflare) ; 


m=(tflare/t3)+1; 

a=t_peak+1; 

G=gammaz(a); % Find ’G’ with the help of ‘gammaz’ function 
fa=(t.*(t_peak) .*exp(-t))./G;% gamma equation 

J=(fa.*J_max) *(1/max(fa)); % Normalized to the total Radiant Intensity 
$J_max 

Js=J*ratwtl; 

fel 2 —size(gec)- 

el=Js(l, 2272); 

c=(zeros(1,m),¢cl1]; 


nie yA 

y=l:length(c) ; 

Yee, 
else 

y=[y,zeros(1, (length(c) -length(y)))]+c; 
end 
BEpL=aO7eS:( (rengen(6)-1)*t3) > 
eval([’save data’,num2str(z),’-.mat’]); 
eval([’load data’,num2str(z),’.mat’]); 


plCtGeEplac) 

Grid! on 
xlabel(’‘time (sec) ’) 
ylabel(’J(W/sr) ’) 
hold on 

end 


figure (2) 
Eeomp-07ts  (( lengenty) =) 73) 
DUOEGEComMp: y ) 

Grldor 

xlabel(’time (sec) ’) 
ylabel(’dJ(W/sr)’) 


% Calculation of radiant intensity for the plume 
Wplb=dquadl1(’plankpl’,lambdal,lambda2);%Emittance between the two lambda 
$( W cm*(-2) ) 

a=1.50; minor axis (m) 

b=3.50; %major axis (m) 

Apl=pi*a*b; % (Area) 

e=0.9; t%emissivity 

s=5.67*10%(-12); %Stefan-Boltzmann constant 

Wplume=s*Tp1%*4; %*Stefan-Boltzmann(total emittance over all wavelengths) 
s1=5.67*10%(-8);% Stefan-Boltzmann constant (m) 

ratWpl=Wplb/Wplume; 

Jpl=((e*si*Tpl%*4)/pi)*Apl*ratWpl; *tRadiant Intensity between the two 
$lambda(W/sr) 

Jplume=Jpl.*ones(1,length(y) ); 


¥ Kam KK KKK S/N Ratio KKKKK KK 
Ratio=Jplume./y; 


figure (3) 
plot (tcomp, Ratio) 


a2 


Giecd. on 
xlabel(’ Pyrophoric Function Time (sec) ‘) 
ylabel(’S/N Ratio’) 


figure (4) 
ror x=l;> lengthier); 
tflare=tf1l (x); 
if length(t_peak1)==1; 
t_peak==t_peakl1; 
else 
t_peak=t_peakl1 (x) ; 
end 
Peelengen lo) =— 1 
(espe els 
else 
Rot— hp 
end 
if length(rmaxp) ==1; 
rmax==rmaxp; 
else 
rmax=rmaxp (x) ; 
end 


% Calculation of radiant intensity for the pyrophoric 

e=-0.9; %emissivity of Pyrophoric 

s=5.67*10%(-12);% Stefan-Boltzmann constant (cm) 

Wpfb=dquad2 (’plankpf’, lambdal,lambda2) ;%Emittance between the two lambda 
%$( W cm*(-2) ) 

Wflare=s*Tpf*4; %Stefan-Boltzmann(total emittance over all wavelengths) 
Afl=4*pi*(rmax)*2; % Area Pyrophoric at t_peak 

s1=5.67*10%(-8);% Stefan-Boltzmann constant (m) 
gemax=((e*sl*Tpt-4) /(4*pi) )*AEl; $$ Total Radiant. tneensityiy/7sr) 
ratWfl=(Wpfb/Wflare); 


m=(tflare/t3) +1; 

a=t_peak+1; 

G=gammaz(a); % Find ‘G’ with the help of ‘gammaz’ function 
fa=(t.*(t_peak) .*exp(-t))./G;% gamma equation 
J=(fa.*J_max) *(1/max(fa)); % Normalized to the total Radiant Intensity 
$J_max 

Js=J*ratwfl; 

[r1,r2]=size(Js); 

el=as(1,2:r2); 

c2=[{zeros(1,m),cl1]; 

Jplume=Jpl.*ones(1,length(c2)); 

c=Jplume./c2; 

Boletos = ( {length (c)=1)*t3) ; 


eval([’save datas’,num2str(x),’.mat’]); 
eval(l’ load datas*,num2str(x),”’.mat‘’ |): 
plortepl, c) 

Grid on 


xlabel(’time (sec) ’) 

ylabel(’S/N Ratio’) 

Hote on 

end 

EEEEEEEEEEEEEEEEESEEEE SECONG Part 

else 

Npix=input(’Set the number of pixels (it must be an integer value with 
half an odd number): '); 

tspf=input(’Set the time (sec) that you want to create the images: '); 


Se 


Tpl=input(’Set the Temperature of plume (Kelvin): '); 

Tpf=input(’Set the Temperature of Pyrophoric (Kelvin): ’); 
lambdal=input(’Set the value of the lambdal (um): '’); 
lambda2=input(’Set the value of the lambda2 (um): '); 

tl=input(’Set the initial value of the time for the Pyrophoric (sec): 
oN 

t2=input (‘Set the final value of the time for the Pyrophoric(sec): '); 
t_peak=input(’Set the value of the t_peak (sec): '); 

rmax=input(’Set the value of the rmax (m): '); 

t3=input(’Set the increment step between tl and t2: '); 


% Calculation of radiant intensity for the plume 
Wplb=dquadl (’plankpl’,lambdal, lambda2); 

a=1.5; minor axis (m) 

b=325-" “smager axis ~(m) 

Apl=pi*a*b; % (Area) 

e=0.9; %emissivity 

s=5.67*10%(-12); %*Stefan-Boltzmann constant 
Wplume=s*Tpl1%*4; %total emittance over all wavelengths 
s1=5.67*10%(-8);% Stefan-Boltzmann constant (m) 
Jplb=((e*s1l*Tpl*4)/pi) *Apl* (Wplb/Wplume); %Radiant Intensity between the 
two lambda (W/sr) 


BOV=((pi*2) 7/150); 
Range=1.8* (10%3); %(m) 
Rfield=FOV*Range; %(m) 


%I create the image matrix for the Plume 
Cpl=zeros (Npix+1,Npix+1); 
forest l=—Norx/ 2 +: Nb x7 2; 
EOt ple —-Noix/ 2: Note 2, 
1pix1=11+((Npix/2)+1); 
7pixi=7 1+ ( (Npix/ 2) +1); 
l= (ie ee (ay Reever). 2) 471). O27 (O/REvela) 2)? 
if wb<=(Noix+ 1) 72; 
CplGipixl, jpixl) =Jpib:% I assign the value Jplb to Ehe pixels 
inside the matrix 
else 
Cpl(ipixl1,jpix1)=0;% I assign the value 0 to the pixels outside 
the matrix 


end 
end 
end 
suml=0; 
fOmet— 1: (Noa 1) 
for j=l. (Npixt)) : 
suml= suml+ Cpl(i,j); 
end 
end 
newCpl= Cpl/suml *Jplb; 
Cpln=newCpl ; 


% Calculation of radiant intensity for the pyrophorice 

e=0.9; %emissivity of Pyrophoric 

s=5.67*10%(-12);% Stefan-Boltzmann constant (cm) 

Wpfb=dquad2 (‘plankpf’,lambdal, lambda2);%Emittance between the two lambda 
%$( W cm*(-2) ) 

Wflare=s*Tpf*4; %total emittance over all wavelengths 

Af1l=4*pi*(rmax)*2; % Area Pyrophoric at t_peak 

s1=5.67*10%(-8) ;% Stefan-Boltzmann constant (m) 
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J_max=((e*sl*Tp£“4)/(4*pi))*Af1l; % Total Radiant Intensity(W/sr) 
ratwWwfl=(Wpfb/Wflare) ; 


(ei eM Oeil te Bie 
vel=rmax./t_peak; 
BaMax estat. ~Vele 
rnorm=r_max_t/rmax; 


% I create the time distribution(gamma distribution) for a=t_peak+1 
a=t_peak+1; 

G=gammaz(a); % Find ’G’ with the help of ‘gammaz’ function 
fa=(t.*(t_peak) .*exp(-t))./G;% gamma equation 
J=(fa.*J_max) * (1/max(fa)); %* Normalized to the total Radiant Intensity 
$J_max 


% I create the normalized space distribution 


s=1; % standard deviation 

fOr exX=C leek ots; 

fad= (exp (-(rnorm- (tx/t_peak) ) .*2./(2*s%2)) ./(sqrt(2*pi)*s)); %Gaussian 
$distribution 

1 Se Ks 


fa2=(12%* (t_peak) *exp(-12)) /gammaz(t_peak+1); 
J2=(fa2*J_max)*(1/max(fa)); % Normalized to the total Radiant 
$Intensity J_max 

02=(0.5) .*erfc( (-tx/ (t_peak*s)) ./sqrt(2)); % Integral from zero to 
infinity for a Gaussian distribution with mean at tx/t_peak 
scrJl=fad*J2/Q2; %*The weighted values for every Gaussian distribution 
$with mean at tx/t_peak 

end 


ber 2-1: beneeth(tspe) ; 

ts=tspf (z); 

havy= (exp (—(rnorm—(lts/t peak) )).°2./ (2*s*2)) ./ (sqrt (2*pi) *sjoe 
fs— Cs; 

fas=(ls*(t_peak) *exp(-ls) ) /gammaz(t_peak+1) ; 
wis=(fas*J max) *(1/max(fa)): 

@©s—( 075). *erfc((-—ts/ (t_peak*s))./sqrt(2));: 

scrJ3=fay*Js/Qs; 


$I create the integrated Pyrophoric model 

Np=30; 

fie ei —size(rnorm ; 

romax=max (rnorm) ; 

Nro=cl1-1; 

Intscr2=1:Nro+l; 

ro=1:Nro+1; % ro is the distance from the center of projection 
Dro=romax/Nro; 

for J=1:Nro+l1; 


sum=0; 
ro(J)=(J-1).*Dro; % for every J pick up one ro 
xmax=sqrt ((romax) .*%2-(ro(J)).%2); 


Dx=xmax./Np; 
for n=1:Np+i1; 
 x=(n-1) .*Dx; 
moan (sc. 24 (Pou) ) 22) > 
ime tax (rnorm) ) Cl: 
k=round(k1) ; % index k goes to closest integer 
te k> Cl 
a=0; 
sum=sum+a; 


re 


elseif k==0; 


a=Q; 
sum=sumt+a; 
else 
sum=sum+scrJd3(k); 
end 
end 
inmeser2 (0 )=Ssum- 
end 


%I create the image matrix for the Pyrophoric model 
POV=t pr 2) / Leo ye 
Range=1.8*(10%3); %(m) 
Rfield=FOV*Range; %(m) 
[ri,cl}=size(rnorm); % cl takes the value of size(r_max_t) 
Cpf=zeros (Npix+1,Npix+1); 
£Om 1=-Npis7 2: 1: Npix/ 2; 
fOr J=—-Npix/ 2:1: Npix/ 2; 
ipax=1+( (Npix/2)41); 
Jpix=j+( (Npix/2)+1); 
Ne=Sorte(i. 24+3.°2) 
R=(Rfield/ (Npix+1))*Nr; 
L1=(R/ (max (rnorm) *rmax) ) *cl1; 
L=round(L1); 
reed 
Cpf(ipix,jpix) =0; 
elseif L==0; 
Cpt (ipix, jpix) =0; 
else 
Spl (lpi, jeio) =Lntserz (Lb); 
end 
end 
end 
Jts=Js*ratwE£l; 
sum=0; 
for  t—1 = (Nor) 
for j3=1: (Npix+1); 
Sums sume Corti, 3): 
end 
end 
newCpf= Cpf/sum *Jts; 
Cpfn=newCpf; 


% keke ke keke KE KK S/N Ratio kak KkKE KK 
Ratio=Jplb/Jts 


eval({’save data’,num2str(z),’.mat’}); 
end 


for 2— length (tspt) : 


eval({’load data’ ,num2str(z),’.mat’]); 
figure (z) 

image (Cpfn, ‘CDataMapping’, ‘scaled’ ) 
end 


figure (z+1) 
image (Cpln, ‘CDataMapping’, ‘scaled’ ) 


end 
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Program 2 Plankpf.m 


% This function file calculates the Spectral Radiant Emittance (SRE) 
% for the pyrophoric at temperature Tpf(Kelvin). 


function Y = plankpf£ (Lambda, T) 

giliebal, Tipe 

c1=3.7415e4; 

c2=1.4388e4; 

VWe=-F Cle, (Gambda ~. 5 2 ae > (CZaeeatiambda = * = Tpt)) —1) ); 
% Y will have units of watts/(cm*2-micron) 


Program 3 Plankpl.m 


% This function file calculates the Spectral Radiant Emittance (SRE) 
% for the plume at temperature Tpl(Kelvin). 


function Y = plankpl (Lambda,T) 

global Tpl 

c1=3.7415e4; 

c2=1.4388e4; 

We Clee, (Lambda .“S 7* (exp(c2 ./(Lambda .<* Tpl))-1))- 
% Y will have units of watts/(cm*%2-micron) 


Program 4 Dquadi.m 
% DQUAD1 Numerically integrates an expression using QUAD8. 
function [1I]=dquadl (Fun,al,a2,a3,a4) 


var = 'x’; 

tol=le-3; 

ie nargin<3 , 

help dquadl 

return 

end 

Lid “igyekelephel ae Si 
xminz=al; 
xmax=a2; 

elseif nargin == 4, 
tee (Si Ste tata: a jue 
var=al; 
xmin=a2; 
xmax=a3; 
else 
San— ele 
xmax=a2; 
tol=a3; 
end 

else 
var=al; 
Sania + 
xmax=a3; 
tol=a4; 

end 

IT=quad8 (’plankpl’,xmin,xmax,tol); 


a7 


Program 5 Dquad2.m 
% DQUAD2 Numerically integrates an expression using QUAD8. 
function [1] =dquad2 (Fun,al,a2,a3,a4) 


Vatuen mex 
tol=le-3; 
ie Naboin< 3, 
help dquad2 
return 
end 
if nargin == 3, 
smin=al ; 
xmax=a2; 
elseif nargin == 4, 
ie (1SStr(al)), 
var=al; 
SM maiz 5 
xmax=a3; 
else 
xmin=al; 
xmax=a2; 
tol=aa- 
end 
else 
var=al; 
xmin=a2; 
xmax=a3; 
fol=a4 
end 
I=quad8 (‘plankpf’,xmin,xmax,tol); 


Program 6 Gammaz.m 


% GAMMAZ gamma function that also allows for complex arguments, unlike 


% MATLAB’s GAMMA function. For real arguments, GAMMAZ gives the 
% same results as GAMMA (within numerical error). 

Funce1om [£] = gammaz (Zz); 

[we rm], = size(z)- 

£ = Zeros (m,n) > 


cof = [76.18009172947146 
=Goascs 22052041677 
24.01409824083091 
=P oy Soe Ze oD 
0.1208650973866179e-2 


~0.5395239384953e-5]; % 6 coefficients in series expansion 

hom ecole —- ism, % do one column at a time 

77 ame (ora i x(@ lS) 

Zp = onesi6, 1) “zz° + [1 lengtm(cof) |) *ones( 1 length(zz))> tevectoerize 

Ser = (cor *(1./ zp) +1. 00000000019 0015) -- 

EMD? = 2245.5 = ((274.-5) . “1004 2225-5)-¢ 

lngamma = -tmp + log(2.5066282746310005*ser./zz); 

£(:,icol) = exp(lngamma) ; 


end 
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APPENDIX C. MATHEMATICAL DERIVATIONS 


1. Solution for t.=a-1 


The pdf of a standard Gamma RV is given by equation (III- 


7) in Chapter III. For convenience it is reproduced below 


l ey, 
foo om 


P(a) 





f(t.a)= 


t>0 
a>0O 


Now by setting 


d 

— f(t,a)=0 

5 f(a) 

with T(a)=(a-1!#0 for a>0 


we get 





| a-2 4-1 4 pe-l f_y\o|=90=—5 
Tay Ma e (-1) 0 


=> (Ge +77 (-le"|=0= a=t+l 


Now, for t=t, 


t =a-l 


oo 


(C-1) 


O=3) 


(C~4) 


Ze Positive side integration of Gaussian distribution 
The pdf of a Gaussian RV is given by equation (III-10) in 


Chapter III. For convenience it iS reproduced below 


l saa : Re eed 
ed a era ao —o<iU<co mean (C-5) 
O<o standard deviation 


The above equation for w=0 and o=1 becomes 


eh —co << X < 00 





(C-6) 


1 
Ol) = 
hoi a pa 


To find the area under the curve from x to infinity the 


following integral must be evaluated: 


"eV dy (C-7) 


In terms of the complementary error function 





ame = 


O(x) = Ya ere he (c-8) 


Now, if the continuous random variable X has mean wu and 


Standard deviation o then 


OO 


z= (c-9) 





is a Gaussian random variable with w=0 and o=1, and the area 


under the curve from z to infinity 1s 


O(z;0,1) =a fe yh ay (C -10) 


From equation (C-8), we have 


at2)= ert 5) (C-11) 


The area under the curve from zero to infinity is obtained 


by setting x=0 in the previous equation to obtain [5] 


at2)=Yert{ Y. 5) (C -12) 


O01 
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